# This script runs all the analyses from the main paper, 
# with the exception of the heckman probit model,
# which is done in Stata.

rm(list=ls())

# Load libraries
library(readstata13)
library(lme4)
library(merTools)
library(ggplot2)
library(stargazer)
library(gridExtra)
library(grid)

# Load data
dta_res <- read.dta13("1.1_shallow_commitments_data_main.dta")

##################
# Main Models and Table 1
##################

# Define list of inputs for models
inputs <- c()

# Model 1: Using the full set of obs, 1982-2015, dem and conc_obs_lag, no interaction
inputs[1] <- "res_withdraw ~ dem + conc_obs_lag + total + res_events + 
                v2x_gender_diff_rat_lag + t_with + t_with2 + 
                t_with3 + (1 | ccode)"

# Model 2: Keep v2x add SOLS, no interaction
inputs[2] <- "res_withdraw ~ dem + conc_obs_lag + total + res_events + v2x_gender_diff_rat_lag +
                sols_lag + t_with + t_with2 + 
                t_with3 + (1 | ccode)"

# Model 3: Sub v2x with wrindex, no interaction
inputs[3] <- "res_withdraw ~ dem + conc_obs_lag + total + res_events + wri_diff_rat_lag +
                sols_lag + t_with + t_with2 + 
                t_with3 + (1 | ccode)"

# Model 4: Add the interaction to mod1
inputs[4] <- "res_withdraw ~ dem + conc_obs_lag + conc_dem + total + res_events + 
                v2x_gender_diff_rat_lag + t_with + t_with2 + 
                t_with3 + (1 | ccode)"

# Model 5: Add the interaction to mod2
inputs[5] <- "res_withdraw ~ dem + conc_obs_lag + conc_dem + total + res_events + 
                v2x_gender_diff_rat_lag + sols_lag +
                t_with + t_with2 + 
                t_with3 + (1 | ccode)"

# Model 6: Add the interaction to mod3
inputs[6] <- "res_withdraw ~ dem + conc_obs_lag + conc_dem + total + res_events + 
                wri_diff_rat_lag +
                sols_lag + t_with + t_with2 + 
                t_with3 + (1 | ccode)"

# Run models
mod1 <- glmer(inputs[1], data = dta_res, family = binomial(link = "logit"), nAGQ = 0, 
              control = glmerControl(optimizer = "nloptwrap"))
# summary(mod1)

mod2 <- glmer(inputs[2], data = dta_res, family = binomial(link = "logit"), nAGQ = 0, 
              control = glmerControl(optimizer = "nloptwrap"))
# summary(mod2)

mod3 <- glmer(inputs[3], data = dta_res, family = binomial(link = "logit"), nAGQ = 0, 
              control = glmerControl(optimizer = "nloptwrap"))
# summary(mod3)

mod4 <- glmer(inputs[4], data = dta_res, family = binomial(link = "logit"), nAGQ = 0, 
              control = glmerControl(optimizer = "nloptwrap"))
# summary(mod4)

mod5 <- glmer(inputs[5], data = dta_res, family = binomial(link = "logit"), nAGQ = 0, 
              control = glmerControl(optimizer = "nloptwrap"))
# summary(mod5)

mod6 <- glmer(inputs[6], data = dta_res, family = binomial(link = "logit"), nAGQ = 0, 
              control = glmerControl(optimizer = "nloptwrap"))
# summary(mod6)


# Build Table 1

stargazer(mod1, mod2, mod3, mod4, mod5, mod6,
          out = "table1_030120b.tex",
          dep.var.labels = "Withdrawing at least one reservation",
          no.space = TRUE,
          keep = c("dem", "conc_dem", "conc_obs_lag", "v2x_gender_diff_rat_lag",
                   "wri_diff_rat_lag", "sols_lag", "res_events", "total",
                   "t_with", "t_with2", "t_with3", "Constant"),
          order = c("dem",  
                    "conc_dem", 
                    "conc_obs_lag",
                    "total",
                    "res_events",
                    "v2x_gender_diff_rat_lag",
                    "wri_diff_rat_lag", 
                    "sols_lag", 
                    "t_with", "t_with2", "t_with3", "Constant"),
          keep.stat = c("n", "aic"),
          covariate.labels = c("Democracy (t-1)", 
                               "Democracy (t-1) * Report (t-1)",
                               "Report (t-1)",
                               "Initial Reservations", 
                               "Global Withdrawals (t-1)",
                               "Women's Rights (t-1) (V-Dem)", 
                               "Women's Rights (t-1) (CIRI)",
                               "SOLS Change (t-1)",
                               "Time", "Time2", "Time3", "Constant"),
          add.lines = list(c("Years", "1982-2015", "1982-2009", "1982-2005",
                             "1982-2015", "1982-2009", "1982-2005"),
                           c("Countries", "45", "43", "41",
                             "45", "43", "41"),
                           c("Withdrawals", "51", "39", "31", "51",
                             "39", "31")),
          notes.label = "",
          notes = "Standard errors in parentheses",
          font.size = "footnotesize",
          column.sep.width = "0pt")


##############
# Substantive Effect of Democracy (Figure 1)
##############

# Set aside data used in Model 1
dat_mod1 <- dta_res[, c("res_withdraw", "dem", "conc_obs_lag", "res_events", "total",
                        "v2x_gender_diff_rat_lag",
                        "t_with", "t_with2", "t_with3", "ccode", "year")]
dat_mod1 <- na.omit(dat_mod1)

demvals <- c(0, 1)

# Create empty matrix to hold average predicted probabilities 
preds_dem <- data.frame(matrix(NA, 2, 4))
colnames(preds_dem) <- c("dem", "mean", "upr", "lwr")
preds_dem$dem <- c(0, 1)

# Create empty matrix to hold simulated predicted probabilities
all_sims_dem <- data.frame(matrix(NA, 10000*2, 2))
colnames(all_sims_dem) <- c("dem", "sims")
all_sims_dem$dem <- rep(c(0,1), each = 10000)

# Loop through values of dem, 
# save simulated predictions in matrix sims
for (d in demvals) {
    
    # alter values of dem, conc_obs_lag, conc_dem
    dat_mod1$dem <- d

    # compute predicted probabilities for all obs
    tmppreds <- predictInterval(mod1, newdata = dat_mod1, seed = 1234,
                                type = "probability", level = .95,
                                n.sims = 10000, returnSims = TRUE,
                                stat = "mean")
    
    # extract simulated predicted probabilities
    sims <- attributes(tmppreds)
    sims <- as.data.frame(sims[-c(1:3)])
    sims <- exp(sims)/(1+exp(sims)) # convert to probability scale
    sims_pred <- as.data.frame(colMeans(sims))
    all_sims_dem$sims[all_sims_dem$dem == d] <- sims_pred$`colMeans(sims)`
    
    # save predicted probabilities
    preds_dem$mean[preds_dem$dem == d] <- mean(sims_pred$`colMeans(sims)`)
    preds_dem$upr[preds_dem$dem == d] <- quantile(sims_pred$`colMeans(sims)`, 0.025)
    preds_dem$lwr[preds_dem$dem == d] <- quantile(sims_pred$`colMeans(sims)`, 0.975)

  }

# Run t.test to verify differences between dem/nondem
t.test(all_sims_dem$sims ~ all_sims_dem$dem) # p < .01

# Figure 1: Graph histograms for dem sims
subeffect_dem <- ggplot(data = all_sims_dem, aes(x = sims)) + 
  geom_histogram(data=all_sims_dem[all_sims_dem$dem == 0,], 
                 aes(fill = "Nondemocracy"), # color = 'grey', fill = 'grey',
                 bins = 100) +
  geom_histogram(data=all_sims_dem[all_sims_dem$dem == 1,],
                 aes(fill = "Democracy"), # color = 'black', fill = 'black',
                 bins = 100) +
  labs(x = "Probability of Withdrawal",
       title = "Substantive Effect of Democracy",
       y = "Frequency") +
  scale_fill_manual(name = "", values = c("black", "grey")) + 
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 10),
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 8),
        text=element_text(family="serif"),
        legend.position=c(0.85,0.8))


# Save Figure 1
ggsave(
  filename = "subeffect_dem.pdf",
  plot = subeffect_dem,
  device = "pdf",
  width = 5.14,
  height = 3,
  dpi = 600
  )


##############
# Substantive Effect of Report Consideration (Figure 2)
##############

# Set aside data used in Model 4
dat_mod4 <- dta_res[, c("res_withdraw", "dem", "conc_obs_lag", "conc_dem",
                        "res_events", "total", "v2x_gender_diff_rat_lag", 
                        "t_with", "t_with2", "t_with3", "ccode", "year")]
dat_mod4 <- na.omit(dat_mod4) # 860 obs

demvals <- c(0, 1)
cvals <- c(0, 1)

# Create empty matrix to hold predicted probabilities 
# (predictions for each of the four combinations of dem and conc_obs_lag)
preds_dc <- data.frame(matrix(NA, 4, 5))
colnames(preds_dc) <- c("dem", "conc_obs_lag", "mean", "upr", "lwr")
preds_dc$dem <- c(0, 0, 1, 1)
preds_dc$conc_obs_lag <- c(0, 1, 0, 1)

# Create empty matrix to hold simulations 
all_sims_co <- data.frame(matrix(NA, 10000*4, 3))
colnames(all_sims_co) <- c("sims", "dem", "conc_obs_lag")
all_sims_co$dem <- rep(c(0,1), each=20000)
all_sims_co$conc_obs_lag <- rep(rep(c(0, 1), each=10000), 2)

# Loop through values of conc_obs_lag and dem, 
# save simulated predictions in matrix sims
for (c in cvals) {
    for (d in demvals) {
      
      # alter values of dem, conc_obs_lag, conc_dem
      dat_mod4$dem <- d
      dat_mod4$conc_obs_lag <- c
      dat_mod4$interaction <- d*c
      
      # compute predicted probabilities for all obs
      tmppreds <- predictInterval(mod4, newdata = dat_mod4, seed = 1234,
                                  type = "probability", level = .95,
                                  n.sims = 10000, returnSims = TRUE,
                                  stat = "mean")
      
      # extract simulated predicted probabilities
      sims <- attributes(tmppreds)
      sims <- as.data.frame(sims[-c(1:3)])
      sims <- exp(sims)/(1+exp(sims)) # convert to probability scale
      sims_pred <- as.data.frame(colMeans(sims))
      all_sims_co$sims[all_sims_co$dem == d & 
                         all_sims_co$conc_obs_lag == c] <- sims_pred$`colMeans(sims)`
      
      # save predicted probabilities
      preds_dc$mean[preds_dc$dem == d & 
                      preds_dc$conc_obs_lag == c] <- mean(sims_pred$`colMeans(sims)`)
      preds_dc$upr[preds_dc$dem == d & 
                     preds_dc$conc_obs_lag == c] <- quantile(sims_pred$`colMeans(sims)`, 0.025)
      preds_dc$lwr[preds_dc$dem == d & 
                     preds_dc$conc_obs_lag == c] <- quantile(sims_pred$`colMeans(sims)`, 0.975)
    }
}

# T-tests for nondemocracies
t.test(sims ~ conc_obs_lag, data = all_sims_co[all_sims_co$dem == 0,])

# Calculate percent increase
new_val <- preds_dc$mean[preds_dc$dem == 0 & preds_dc$conc_obs_lag == 1] 
old_val <- preds_dc$mean[preds_dc$dem == 0 & preds_dc$conc_obs_lag == 0] 
(new_val - old_val)/old_val

# Figure 2: Create histograms
se_concobs_d0 <- ggplot(data = all_sims_co, aes(x = sims)) + 
  geom_histogram(data=all_sims_co[all_sims_co$dem == 0 & all_sims_co$conc_obs_lag == 0,],
                 aes(fill = "No Report"),
                 bins = 100) +
  geom_histogram(data=all_sims_co[all_sims_co$dem == 0 & all_sims_co$conc_obs_lag == 1,], 
                 aes(fill = "Report"),
                 bins = 100) +
  labs(x = "Probability of Withdrawal",
       title = "Nondemocratic States",
       y = "Frequency") +
  xlim(0,0.5) +
  scale_fill_manual(name = "", values = c("grey", "black")) + 
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 10),
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 8),
        text=element_text(family="serif"),
        legend.position=c(0.7,0.8))

se_concobs_d1 <- ggplot(data = all_sims_co, aes(x = sims)) + 
  geom_histogram(data=all_sims_co[all_sims_co$dem == 1 & all_sims_co$conc_obs_lag == 0,],
                 aes(fill = "No Report"),
                 bins = 100) +
  geom_histogram(data=all_sims_co[all_sims_co$dem == 1 & all_sims_co$conc_obs_lag == 1,], 
                 aes(fill = "Report"),
                 bins = 100) +
  labs(x = "Probability of Withdrawal",
       title = "Democratic States",
       y = "Frequency") +
  xlim(0,0.5) +
  scale_fill_manual(name = "", values = c("grey", "black")) + 
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 10),
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 8),
        text=element_text(family="serif"),
        legend.position=c(0.7,0.8))

subeffect_concobs <- grid.arrange(se_concobs_d0, se_concobs_d1, nrow=1, ncol=2, 
             top = textGrob("Substantive Effect of Report Consideration",
                            gp = gpar(fontsize = 10, fontfamily="serif"), vjust = 0.5,
                            hjust = 0.5))

# Save Figure 1
ggsave(
  filename = "subeffect_concobs.pdf",
  plot = subeffect_concobs,
  device = "pdf",
  width = 6,
  height = 3.5,
  dpi = 600
)
